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Abstract 


A new approach utilizing Newton Method and Homotopy Analysis Method 
(HAM) is proposed for solving nonlinear system of equations. Accelerating 
the rate of convergence of HAM, and obtaining a global quadratic rate of 
convergence are the main purposes of this approach. The numerical results 
demonstrate the efficiency and the performance of proposed approach. The 
comparison with conventional homotopy method, Newton Method and HAM 
shows the great freedom of selecting the initial guess, in this approach. 


Keywords: Homotopy Analysis Method; Zero order deformation equations; 
Control convergence parameter; Newton’s method; Iterative method; Multi- 
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1 Introduction 


Solving algebraic and transcendental equations is an interesting mathemat- 
ical problem that has been occupied an important place in mathematical 
history. This problem arises in different applications of mathematics in sci- 
ences and engineering. Analytical solution of this problem is reserved to 
a small category of equations. For this reason and the exigencies of those 
increasing applications, from the beginning of era of electronic computing 
numerical methods of these problems have been progressed. 
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Actually there is a vast group of conventional methods to solve algebraic 
and transcendental equations, but yet there exist enormous difficulty due to 
local convergence of these methods that make the new research inevitable. 
Particular numerical solution of system of nonlinear equations is realized by 
different methods. A traditional method is Newton method that can have 
quadratic order of convergence, but the convergence is local [16]. There 
is a variety of modified Newton methods which make a global convergence 
possible [16]. Many new one-step and multi-step methods are used to solve 
these system of equations (for more details one can refer to [4, 7, 8]). There 
is also acceleration methods and multi-step methods but these methods are 
also very dependent to initial guess and have local convergence in the most of 
the cases [16]. Recently the homotopy method using the notion of homotopy 
and functional series are applied to solve the system of nonlinear equations 
[1, 3, 6, 15, 11, 17]. Some methods are very suitable, but in practice they 
need to solve a system of differential equations with initial conditions [14]. 
One of the most important of Homotopy methods which is principally used 
for solving the nonlinear differential equations is Homotopy Analysis Method 
(HAM), that can be applied for solving nonlinear equations, but it is normally 
slow with local convergence [14]. In this paper a combination of Newton 
Method and HAM is considered to solve the algebraic and transcendental 
system of equations with the aim of improving the both mentioned methods, 
in view of local convergence and the rate of convergence. The results of 
proposed method will be compared with other methods. 

The organization of the paper is as follows. In Section 2 a concise descrip- 
tion of the Newton Method, the Homotopy Method are presented. In Section 
3 the fundamental of HAM and proposed approach is discussed. In Section 4 
the numerical results for 3 methods are given and compared. Finally Section 
5 ends the paper with conclusion and discussion. 


2 Description of problem and the methods 


Consider the following nonlinear algebraic or transcendental system of equa- 
tions 

F(x) =0, F = (fi, fas-*+sfn); (1) 
where F : D Cc R” > R”, that D is an open region in R” and F € C1(D) 
such that F(x) = 0. The vector X is called the zero of F' or the solution of the 
equation (1). Recalling that the Newton Method for solving (1) is formulated 
as follows 


x(P+1) — x() _ [DF (x) -1 F(x), k=0,1,2,--:, (2) 


where DF is the Jacobian matrix of F and x) is an initial guess of &. 
For more details see [16]. The Newton method is a suitable technique for 
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differentiable functions. In general, the rate of convergence is quadratic in a 
neighborhood of the solution X, with local convergence property. As a second 
choice for solving (1), the homotopy method for the system of nonlinear 
equation is recalled [6]. The Homotopy function 


H: (0,1) xR” > R”, 
is defined by 


H(q,x) = F(x) + (1— g)(F(%) — F(x) (3) 
= F(x) + (@-1) Fe), (4) 


here x) is an initial guess of X and q is called Homotopy parameter or 
embedding parameter. Obviously, at gq = 0 and q = 1, 


H(0,x) = F(x) -— F(x), = H(1,x) = F(x). 


If g increases from 0 to 1 then the function H(q, x) varies continuously from 
F(x) — F(x) to F(x). In topology, such a kind of continuous variation 
is called deformation. The function H respect to parameter qg, provides us 
a family of functions that can lead from the known value x), to solution 
x. The function H is a Homotopy between H(0,x) = F(x) — F(x)) and 
H(1,x) = F(x). Accepting that ¢ : [0,1] > R", x = ¢(q) is a unique solution 
of the equation 

H(q, x) =0, qe (0, 1] ’ (5) 


or 
H(,e(g)=90, ge (0,1). (6) 


The set {¢(q)|O < q < 1} can be viewed as a family of parameterized curves 
respect to q in R” from ¢(0) to ¢(1) = xX. The solution X of F(x) = 0 can be 
obtained by solving the following system of equations 


¢'(9) = -[F(O@)"F((0), = OSG <1, 


with the initial condition (0) = x), where J(¢(q)) is jacobian matrix of H 
respect to x [6]. This method will be reffered as HM. 


3 HAM combined with Newton method 


The Homotopy Analysis Method (HAM) is proposed by Liao [2]. In this 
method one introduces a homotopy function for solving (1). To be more 
precise, the following homotopy function is considered: 


H[q, 4(q)] = (1 - L[o(q) — x] + aN [d(Q)] , (7) 
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where q € [0,1] is an embedding parameter and ¢(q) is a function of g, and 
x) € R” is an initial estimation of X, the solution of (1). Also, N is a 
nonlinear operator and CL is a linear operator and 


N(x) = F(x). (8) 
If ¢ = 0 and q = 1, then considering ¢(0) = x, yields 
Hq, $(9)]|,-9 = £190) — x] =0, (9) 


and 


Ha, 9(@)l|,. =NIOY), oo) =®- 


By using (9), the vector 
ol) =x, 


is obviously the solution of the equation 
Ha. 6(@)]|,-4 =9- 


As the embedding parameter qg increases from 0 to 1, the solution ¢(q) of 
equation 


Hq, o(q)] = 0, 


depends upon the embedding parameter q and varies from initial approxima- 
tion x) to the solution ¥ of equation (9). Now by using homotopy function 
(7) we construct a family of equations 


(1—@)L[e(q) - x] =aNV(0@), = € [0,1], (10) 
subject to the initial condition 


o(0) =x , (11) 


Consider equation (1) and let A be a non-singular matrix which will be 
determined later. We construct following deformation equation that is called 
zeroth-order deformation equation: 


(1 — g)A(o(q) — x) = gF(¢(q)) . (12) 


Suppose X is solution of F(x) = 0 and the sequence {xo exist with the 
iEN 


following property 


and 


A new approach for solving nonlinear system of equations ... 61 


Next a function ¢ : [0,1] — R” is defined as follows 
x= 4(q)= >> xq", = g € (0,1, 


Subject to 


By differentiating (12) with respect to q, the following equation is obtained: 


—A(o(q) - x) + (1 a(A+-0(4)) = F(¢(q)) + <_F(o4a)) (15) 


Putting g = 0 in (15) yields 


d 
AT o(@)| = F(6(0) (16) 
q q=0 
Matrix A being non-singular, it deduces 
d -1 
“| = A-1F(¢(0) (17) 
q a 
On the other hand, 
d co 
eae = mxi™ gm 
aq? 2 q 
Then 
—4(g)} =x = ATF(x) (18) 
dq 526 


The equation of (15) is called first-order deformation equation. By differen- 
tiating equation (15) with respect to q, the following equation is obtained 


qq?) ( q) dq o(q) 
d d? 


= 2g TP) +993 


F(9(q)) - (19) 
Putting g = 0, the second-order deformation equation is obtained as follows 
—2Ax") 4+ 2Ax® = 2D, F(x)x® , (20) 


or 
x?) = (AD, F(x) + 1D)x® . 
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By repeating the same procedure the m-th order deformation equation can 
be obtained. Indeed, the following proposition can be proved. 


Proposition 3.1. Jf F : R" > R”, and F ¢ C™(R"), A € R"™” a given 
matriz, and 


z=¢d(g)= > a", =o: [0,1] 9”, 
m=0 


(1 — q)A(o(q) — &) = gF(9(q)). 


where @ is an analytic function, then 


ii ee 
Alig —y gD) — ——— —_ ee 21 
(a Nine ) (m—1)! dg (0(q)) hoe (21) 
where 
_ so m<l1 
Ae Ny 0.W. 
Ifm > 2 and A be a nonsingular matrix then 
1 ia 
(m) _ g(m-l) | Av} F 22 
x x Gut)! dg”! (0(q, 2)) sag (22) 


The equation (22) is called m-th order deformation equation. 


For solving system of algebraic equations in general one can use the above 
equations to determine the vectorial terms x of ® = 772)x™, ie. the 
following equations. 


-1 (0) = 
g(r = | AFA) ae gs 
xO") + hy A! dinar F(6(9) Meee 
q=0 


In practice, one can obtain a finite number of x, Then by considering 
partial sum of above series one can determine ¢(1) approximately by a K*” 
order partial sum as follows: 


K=P(1) 2 xO FxO 4.4% , 


Unfortunately, the homotopy series 
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$(q) = S> xq” , 


m=0 


may be divergent at g = 1. To overcome this restriction, Liao [14] introduced 
an auxiliary parameter h # 0 to construct a kind of deformation equations 
based on 


(1 — q)A(o(q,h) — x) = qhF(6(q,h)) , 


where 
CO 


(a, h) = >> x (h)a™, 


m=0 


the vectors xm) are dependent on h. In particular if series is convergent for 
at least one h, it is deduced [9], 


= S> x (A), d(0,h) =x, 6(1,h) =. 


m=0 


Therefore, the equation (23) is transformed to 


4, (0) = 
xn) = J RAMP) es.) 
m—1 
x09 + GaerA agar F(0(4 %)) ee 
q=0 


The parameter h is called convergence control parameter. The convergence 
rate and region of series solution depend on the convergent control parameter. 
This parameter provides a convenient way to adjust and control convergence 
region and rate of convergence of series solution given by the HAM. For find- 
ing a suitable h, some approaches are proposed in [2, 5]. The traditional 
approach gives the possibility of estimation a suitable value of h, by plot- 
ting the h-curves (for more details see [14]). Following [9], we use a more 
systematic approach in this work. Consider 


K 
(a,b) = O(a,h) = D> xOMQ” = xO 4 xMq4xMP 4... +x qX . 


m=0 


The value a1, h) is only a function of h, which is denoted by 


~N 


Pe(h) = (1, h) . 


As proved by Liao in general [14], if the series solution converges, then there 
exists at least an ho such that 
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jim ll F(we(ho)) |= 0, 


where denote ||.|| is Euclidian norm in R". Accordingly, we let 
F(x (ho))|] = main |e (h))I1 (25) 
LERy, 


where Ry», is a valid region that lie on a horizontal segment of the h-curves. 
The (ho) is a vector in R” that can be regarded as an approximation of 
x. So, we can apply ~x(ho) as initial point for Newton method, if Newton 
method converges, the desired approximate solution is found, otherwise, after 
some iterations, the result of Newton method is considered as an initial point 
for a new HAM procedure and so on. 

The proof of convergence is an open problem [14] . The numerical exam- 
ples show that proposed method is more efficient than Newton method. 

The proposed HAM is convergent for many examples but this method 
spends a lot of time during each iteration. For accelerating the convergence 
this method, we suggest the combination of HAM and Newton method. At 
the beginning, a new initial point can be obtained by utilizing the proposed 
method, then the process continues by Newton method with this new initial 
point. If Newton method does not converge to solution after some itera- 
tions, the HAM method can be applied again by using this new initial point. 
If DF (x) is non-singular, this matrix is practically profitable as a good 
selection of A, so 

A= DF(x), 


Using the above choice it is observed when h = —1 the first step of the 
homotopy consists of the first iteration of Newton method, in fact, one has 


Ee 1 
x = d(1) = d(1) = bs xa = x 4 4Q) . 
q=1 


where by using (24) 
x = —DF(x)-1 FRM) . 


This result demonstrate the validity of choosing A = DF (x). Application 
and implementation of this hybrid method allow us improving local conver- 
gence of newton method , and choosing x) arbitrary. 


4 Numerical experiments 


In this section, several examples are considered and the numerical results 
for mentioned methods: Homotopy Method(HM), HAM, Newton method 
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Table 1: Numerical results for Example 4.1 with xe) = 1 


Method NI F«O))|| CPU time result 

NHAM 4 4.335133e — 008 3.333667e + 000 Convergent 
Newton 4 1.691234e — 008 7.639678e — 002 Convergent 
HAM 3 5.775537e — 008 2.12925le — 001 Convergent 
HM = 7.457211e — 005 4.041735e — 001 Convergent 


log||F (XN) || 


Number of Iteration 


Figure 1: The graph of In(||F(X)||) for Example 4.1 with 2) =1 


and Newton-HAM (NHAM) are reported. We utilize MATLAB 8. In Tables 
and Figures, the number of iterations (NI), the Euclidean norm of residual 
of government equation and CPU time, are presented. 


Example 4.1. Consider the following equation: 
f(x) = xe” —1=0, (26) 


The function f has at least one zero between 0 and 1. For 2 = 1, the 
numerical results are shown in Table 1. For this initial point all methods are 
convergent, but the Newton method is apparently faster than other methods. 
For 2°) = 10, the numerical results are shown in Table 2. 

In this case, HM method is divergent, Newton method is faster than 
NHAM and HAM and results are more accurate than others. The number of 
iterations for NHAM is less than the others. For 2) = —400, the numerical 
results are shown in the Table 3. In this example NHAM method is conver- 
gent and other methods are divergent. 
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Table 2: Numerical results for Example 4.1 with 2) = 10 


Method NI LF (xO")))| CPU time result 
NHAM 6 2.160101e — 010 3.148814e + 000 Convergent 
Newton 16 5.107026e — 014 2.622956e — 001 Convergent 
HAM 31 3.594237e — 008 2.876829e + 000 Convergent 
HM — 5.790573e + 003 4.081270e — 001 Divergent 


—+— NEWTON} 


\ 1 
5 10 15 20 25 30 35 
Number of Iteration 


Figure 2: The graph of In(||F(X)||) for Example 4.1 with «© = 10 


Table 3: Numerical results for Example 4.1 with 2) = —400 


Method NI F(xO"))]| CPU time result 
NHAM 212 2.403109e — 007 6.452634e + 000 Convergent 
Newton 3 Nan = Divergent 
HAM 100 Infinity - Divergent 
HM - Nan - Divergent 
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log||F(XN)|| 


1 if L 
0 50 100 150 200 250 
Number of Iteration 


Figure 3: The graph of In(||F(X)||) for Example 4.1 with «© = —400 


Table 4: Numerical results for Example 4.2 


Method NI EF (xO"))]| CPU time result 
NHAM 6 2.085579e — 008 2.116125e + 000 Convergent 
Newton 1 NaN - Divergent 
HAM 8 2.247981e — 010 3.105904e + 000 Convergent 
HM — 1.967763e + 009 7.584255e — 001 Divergent 


Example 4.2. Consider following equations: 


fila, y, z,d) = xyz +d—31=0, 
fo(a,y,2,€@) =a+y+24+d-11=0, 
fs(a, y, 2, d) = 2a + 8y +42 + d— 35 =0, 
fa(z,y, 2,4) =a+z2-—y+d-1=0, 


(27) 


where F = | fi fo fs ie 


We know X, = (2,3,5,1) and X, = (2, igo 2) are two solutions of 


F(X) =0. For X = (1,1,1,1) , numerical results are shown in Table 4. 


Newton Method is divergent because det(DF(X)) = 0. But HAM and 
NHAM methods converge, and NHAM is faster than HAM. 


Example 4.3. Consider the following equations: 
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20b —a— NHAM 
—+— NEWTON 
—+— HAM 
—s— HM 


log||F(XN)|| 


L L L L L L 
1 2 3 4 5 6 7 8 
Number of Iteration 


Figure 4: The graph of In(||F(X)||) for Example 4.2 


Table 5: Numerical results for Example 4.3 


Method NI FxO"))]| CPU time result 
NHAM 8 6.567317e — 010 8.196840e + 001 Convergent 
Newton 101 5.030214e + 003 6.042236e + 000 Divergent 
HAM 18 5.830347e — 008 3.106638e + 003 Convergent 
HM - 5.242329e + 002 2.483993e + 000 Divergent 


fii, a,°°° stn) = (3 - 5X )Ly —2%2+1 =0 5 


. sce = (3—1,7.)¢; ‘ : re, ee ; 
a ’ ’ son) = tJXd i a+ — ’ ’ 
fi(x1, £2 En) = (3 — 50;)a;— 2-1 — 21 t1=0, l<i<n 


fn(@1,2,°°+ 5 2n) = (8 - alain — 2%n-1+1=0, 
(28) 

that F = [fi fo ful” 

For n = 50 and X) = (100, 100,--- ,100), numerical results are shown 
in Table 5. 

In this example HAM and NHAM are convergent to the exact solution 
X= (1,...,1), but Newton method is divergent. Also NHAM is faster than 
HAM. Results are shown in Figure 5. 


Example 4.4. Consider the following equations: 


fe(@1, €2,°°* , Un) = 100002, 2441 —1=0, mod(k,2) = 1, 


fe(@1,%2,°++ , Un) = exp(—ap-1) + exp(—x,) — 1.0001 = 0, mod(k, 2) = 0, 
(29) 
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—#— NHAM 


ma NN = | 


log||F(XN)|| 


50 7 
Number of Iteration 


Figure 5: The graph of In(||F(X)||) for Example 4.3 


Table 6: Numerical results for Example 4.4 


Method NI LF (xl"))]| CPU time result 
NHAM 13 1.059758e — 009 5.152483e + 001 Convergent 
Newton 12 1.112530e — 010 1.972995e + 001 Convergent 
HAM 23 2.728830e + 056 1.128194e + 002 Divergent 
HM = 4.211734e — 001 1.092275e + 002 Convergent 
T 
that F = [fa fe sity’ Jal 


For n = 100 and X©) = (1,0,1,0,1,0,1,--- ,0), numerical results are 
shown in Table 6. In this example Newton method, NHAM and HM are 
convergent, but HAM is divergent. Results are shown in Figure 6. 


Example 4.5. Consider the following equations: 


fi(ai, v2) = exp(a1) + 2142 —1=0, 
(30) 
fo(a1, x2) = sin(x1£2) +21, +%2—- l= 0, 
T 
that F = [fi fe eat Jal 
For n = 100 and X) = (1,0,1,0,1,0,1,--- ,0), numerical results are 
shown in Table 7. In this example all the methods are convergent. Results 
are shown in Figure 7. 
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log||F(XN)|| 
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Number of Iteration 


Figure 6: The graph of In(||F(X)||) for Example 4.4 


Table 7: Numerical results for Example 4.5 


Method NI LF (xO"))]| CPU time result 

NHAM 4 1.993082e — 010 1.848117e + 000 Convergent 
Newton 4 1.405720e — 012 2.472726e — 001 Convergent 
HAM 4 1.187309e — 008 1.719704e + 000 Convergent 
HM - 5.368713e — 006 1.380355e + 000 Convergent 


log||F(XN)|| 


-30 f : 
1 2 3 4 
Number of Iteration 


Figure 7: The graph of In(||F(X)||) for Example 4.5 


A new approach for solving nonlinear system of equations ... 71 


5 Conclusion 


In this paper, Newton-HAM (NHAM) applying control parameter h are pro- 
posed for solving systems of nonlinear equations. The results for all examples 
are convergent and also NHAM is faster than Homotopy method. The re- 
sults demonstrate that by choosing a suitable h, HAM and NHAM methods 
are convergent. The numerical results show in general that the proposed 
method is effective and efficient and provides highly accurate results in a less 
number of iterations as compared by other methods. The main advantage 
of NHAM is the relative freedom of choosing initial guess. The appropriate 
proof convergence of NHAM can be continuation of the present work. 
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